##Generate data and maps for USDOL Prisoner Reentry Grant (April 2017)
##Load packages
library(tigris)
##
## Attaching package: 'tigris'
## The following object is masked from 'package:graphics':
##
## plot
library(stringr)
library(xlsx)
## Loading required package: rJava
## Loading required package: xlsxjars
library(readxl)
library(leaflet)
library(maptools)
## Loading required package: sp
## Checking rgeos availability: TRUE
library(ggplot2)
library(ggmap)
library(RColorBrewer)
library(rgdal)
## rgdal: version: 1.2-5, (SVN revision 648)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
## Path to GDAL shared files: C:/Users/sschmid/Documents/R/R-3.3.3/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
## Path to PROJ.4 shared files: C:/Users/sschmid/Documents/R/R-3.3.3/library/rgdal/proj
## Linking to sp version: 1.2-4
##Set working directory to Shapefile folder
setwd("~/R Projects/USDOL_ReentryProject/Shapefiles")
##Read in saved shapefiles
####Michigan cities
cities <- readOGR(dsn = "Michigan_Cities",
layer = "tl_2010_26_place00")
## OGR data source with driver: ESRI Shapefile
## Source: "Michigan_Cities", layer: "tl_2010_26_place00"
## with 630 features
## It has 17 fields
## Integer64 fields read as strings: ALAND00 AWATER00
####DPD precincts
precincts <- readOGR(dsn = "DPD Precincts",
layer = "geo_export_7b552edd-0682-452c-b491-7e6a16b1c277")
## OGR data source with driver: ESRI Shapefile
## Source: "DPD Precincts", layer: "geo_export_7b552edd-0682-452c-b491-7e6a16b1c277"
## with 12 features
## It has 12 fields
####Detroit neighborhoods
neighborhoods <- readOGR(dsn = "Detroit Neighborhoods",
layer = "geo_export_9b67ef6c-9a29-4277-87eb-e2d8eafc5186")
## OGR data source with driver: ESRI Shapefile
## Source: "Detroit Neighborhoods", layer: "geo_export_9b67ef6c-9a29-4277-87eb-e2d8eafc5186"
## with 201 features
## It has 1 fields
##Read in shapefiles from tigris package
####Southeast MI census tracts
tracts <- tracts(state="MI", county=c(99, 125, 163), cb=TRUE)
####Michigan counties
counties <- counties(state = "MI", cb = TRUE)
####ZCTA starting with "48"
zipcodes <- zctas(cb=TRUE, starts_with=c("48"))
##Set working directory to Raw Data folder
setwd("~/R Projects/USDOL_ReentryProject/Raw Data")
##Read in raw data files
####ACS 2014 5-year poverty data by census tract
tractfile <- read_excel("Poverty_CensusTracts_492017.xlsx")
####UWSEM school partner names and geocoordinates
schoolfile <- read.csv("UWSEM School Coordinates.csv", header=TRUE)
####DPD 2015 crime incidents
crimefile <- read.csv("Crime_2015.csv", header=TRUE)
##Set working directory to main project folder
setwd("~/R Projects/USDOL_ReentryProject/")
##Clean data and shapefiles
####Keep only SE Michigan counties in counties shapefile
counties <- counties[(counties$NAME == "Macomb" |
counties$NAME == "Oakland" |
counties$NAME == "Wayne"),]
####Remove non-violent crime data
crimefile <- crimefile[which(crimefile$CATEGORY == "AGGRAVATED ASSAULT" |
crimefile$CATEGORY == "HOMICIDE" |
crimefile$CATEGORY =="ROBBERY"),]
####Change neighborhood names in shapefile to match crime data
neighborhoods$name <- toupper(neighborhoods$name)
neighborhoods$name[neighborhoods$name=="BARTON-MCFARLAND"] <-
"BARTON MCFARLANE"
neighborhoods$name[neighborhoods$name=="SOUTHWEST DETROIT"] <-
"SOUTHWEST DETROIT / MEXICANTOWN"
##Merge data to shapefiles
####Merge ACS poverty rate data to census tract shapefile
merged <- geo_join(tracts, tractfile, "AFFGEOID", "AFFGEOID")
##Duplicate census tract shapefile and remove high poverty census tracts
mergedpoverty <- merged
mergedpoverty <- mergedpoverty[!(is.na(mergedpoverty$poverty)),]
mergedpoverty <- mergedpoverty[!(mergedpoverty$poverty > 30),]
####Merge crime data to census tract shapefile
crimetract <- table(crimefile$CENSUSTRACT)
crimetract <- as.data.frame(crimetract)
crimetract$tract <-crimetract$Var1
mergedtract <- geo_join(merged, crimetract, "NAME", "tract")
##Calculate violent crime rate/50000 residents
mergedtract$crimerate <- (mergedtract$Freq/mergedtract$population)*50000
##Duplicate census tract crime file and remove high crime census tracts
mergedcrime <- mergedtract
mergedcrime <- mergedcrime[!(is.na(mergedcrime$crimerate)),]
mergedcrime <- mergedcrime[!(mergedcrime$crimerate > 840),]
##Merge crime data to Detroit neighborhood shapefile
crimeneighborhood <- table(crimefile$NEIGHBORHOOD)
crimeneighborhood <- as.data.frame(crimeneighborhood)
crimeneighborhood$name <-crimeneighborhood$Var1
mergedneighborhood <- geo_join(neighborhoods, crimeneighborhood, "name", "name")
##Merge crime data to DPD precincts shapefile
crimeprecinct <- table(crimefile$PRECINCT)
crimeprecinct <- as.data.frame(crimeprecinct)
crimeprecinct$PRECINCT <-crimeprecinct$Var1
####Clean precinct names in crime data file to match shapefile
crimeprecinct$PRECINCT <- as.character(crimeprecinct$PRECINCT)
crimeprecinct$PRECINCT[crimeprecinct$PRECINCT=="2"] <- "02"
crimeprecinct$PRECINCT[crimeprecinct$PRECINCT=="3"] <- "03"
crimeprecinct$PRECINCT[crimeprecinct$PRECINCT=="4"] <- "04"
crimeprecinct$PRECINCT[crimeprecinct$PRECINCT=="5"] <- "05"
crimeprecinct$PRECINCT[crimeprecinct$PRECINCT=="6"] <- "06"
crimeprecinct$PRECINCT[crimeprecinct$PRECINCT=="7"] <- "07"
crimeprecinct$PRECINCT[crimeprecinct$PRECINCT=="8"] <- "08"
crimeprecinct$PRECINCT[crimeprecinct$PRECINCT=="9"] <- "09"
####Merge files
mergedprecinct <- geo_join(precincts, crimeprecinct, "name", "PRECINCT")
##Add population data to precinct shapefile - source: 2010 Census
####Create data frame with population data
population <- c("54499", "43701", "66398", "49435", "63701", "46342",
"86006", "74348", "56894", "56596", "75052")
precinct <- c("02", "03", "04", "05", "06", "07", "08", "09", "10", "11",
"12")
precinctpop <- cbind(population, precinct)
precinctpop <- as.data.frame(precinctpop)
####Merge files
mergedprecinctcrime <- geo_join(mergedprecinct, precinctpop, "name", "precinct")
####Calculate violent crime rate/50000 residents for each precinct
mergedprecinctcrime$Freq <- as.numeric(mergedprecinctcrime$Freq)
mergedprecinctcrime$population <-
as.numeric(as.character(mergedprecinctcrime$population))
mergedprecinctcrime$crimerate <-
(mergedprecinctcrime$Freq/mergedprecinctcrime$population)*50000
##Create labels and popups
popupschool <-paste0(schoolfile$school)
zipcodelabel <- paste0(zipcodes$ZCTA5CE10)
popup1 <- paste0("Census Tract: ", merged$NAME, "<br>",
"Population: ", merged$population, "<br>",
"Poverty Rate: ", round(merged$poverty, 2), "%")
popupcity <-paste0(cities$NAME00)
precinctlabel <- paste0("DPD Precinct ", precincts$name)
popuplabel <- paste0(neighborhoods$name)
popup2 <- paste0("Census Tract: ", mergedtract$NAME, "<br>",
"Population: ", mergedtract$population, "<br>",
"Violent Crimes: ", mergedtract$Freq, "<br>",
"Violent Crime Rate: ", round(mergedtract$crimerate,0))
popuphood <- paste0("Neighborhood: ", mergedneighborhood$name, "<br>",
"Violent Crimes: ", mergedneighborhood$Freq)
popupprecinct <- paste0("Precinct ", mergedprecinctcrime$name, "<br>",
"Violent Crime Rate: ",
round(mergedprecinctcrime$crimerate, 0))
##Create color palettes for crime and poverty rates
pal <- colorNumeric(palette="Reds",
domain = merged$poverty, na.color="Gray")
qpal <- colorBin(palette="Reds",
domain = mergedtract$crimerate, na.color="Gray",
bins=8)
palhood <- colorBin("Reds", domain=mergedneighborhood$Freq, bins=4)
palprecinct <- colorBin("Reds", domain=mergedprecinctcrime$crimerate)
##Create leaflet map
map1<-leaflet() %>%
addProviderTiles(providers$Stamen.TonerLite) %>%
setView(lng=-83.170268, lat=42.373568, zoom=10) %>%
addPolygons(data=counties, group="Census Tract Poverty Rates",
fillColor="Gray", color="black",
fillOpacity = 0, weight=5) %>%
addPolygons(data = merged, group="Census Tract Poverty Rates",
fillColor = ~pal(poverty),
weight=0.3, fillOpacity = 0.7, color="Gray",
popup = popup1) %>%
addPolygons(data=counties, group="Census Tract Violent Crime Rates",
fillColor="Gray", color="black",
fillOpacity = 0, weight=5) %>%
addPolygons(data = mergedtract, group="Census Tract Violent Crime Rates",
fillColor = ~qpal(crimerate),
weight=0.3, fillOpacity = 0.7, color="Gray",
popup = popup2) %>%
addPolygons(data=counties, group="DPD Precinct Violent Crime Rates",
fillColor="Gray", color="black",
fillOpacity = 0, weight=5) %>%
addPolygons(data = mergedprecinctcrime,
group="DPD Precinct Violent Crime Rates",
fillColor = ~palprecinct(crimerate),
weight=0.3, fillOpacity = 0.7, color="Gray",
popup = popupprecinct) %>%
addLayersControl(baseGroups = c("Census Tract Poverty Rates",
"Census Tract Violent Crime Rates",
"DPD Precinct Violent Crime Rates"),
overlayGroups = c("City Borders",
"Zip Codes",
"DPD Precincts",
"Detroit Neighborhoods",
"High Poverty Census Tracts Only",
"High Crime Census Tracts Only",
"Returning Residents Center",
"UWSEM Schools"),
position="bottomleft",
options = layersControlOptions(collapsed = FALSE)
) %>%
addPolygons(data=cities, group="City Borders", fillColor="Gray",
color="Black", fillOpacity=0, weight=2, dashArray="3",
label=popupcity) %>%
addPolygons(data=zipcodes, group="Zip Codes", fillColor="Gray",
color="Black", fillOpacity=0, weight=2, dashArray="3",
label=zipcodelabel) %>%
addPolygons(data=precincts, group="DPD Precincts", fillColor="Gray",
color="Black", fillOpacity=0, weight=2, dashArray="3",
label=precinctlabel) %>%
addPolygons(data=neighborhoods, group="Detroit Neighborhoods",
fillColor="Gray", color="Black", fillOpacity=0, weight=2,
dashArray="3",
label=popuplabel) %>%
addPolygons(data=mergedpoverty, group="High Poverty Census Tracts Only",
fillColor="Gray", color="Gray", fillOpacity=0.7, weight=1) %>%
addPolygons(data=mergedcrime, group="High Crime Census Tracts Only",
fillColor="Gray", color="Gray", fillOpacity=0.7, weight=1) %>%
addMarkers(lng=-83.170268, lat=42.373568, group="Returning Residents Center",
label="Returning Residents Resource Center") %>%
addMarkers(group="UWSEM Schools", lng=schoolfile$lng, lat=schoolfile$lat,
label=popupschool) %>%
hideGroup("Zip Codes") %>%
hideGroup("City Borders") %>%
hideGroup("DPD Precincts") %>%
hideGroup("Detroit Neighborhoods") %>%
hideGroup("High Poverty Census Tracts Only") %>%
hideGroup("High Crime Census Tracts Only") %>%
hideGroup("UWSEM Schools")
## Warning in qpal(crimerate): Some values were outside the color scale and
## will be treated as NA
##View map
map1